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We numerically study a one-dimensional system of N classical localized planar rotators coupled 
through interactions which decay with distance as l/r a (a > 0). The approach is a first principle one 
(i.e., based on Newton's law) which, through molecular dynamics, yields the probability distribution 
of angular momenta. For a large enough we observe, for longstanding states corresponding to 
TV 2> 1 systems, the expected Maxwellian distribution. But, for a small or comparable to unity, 
we observe instead robust fat-tailed distributions that are quite well fitted with g-Gaussians. These 
distributions extremize, under appropriate simple constraints, the nonadditive entropy S q upon 
which nonextensive statistical mechanics is based. The whole scenario appears to be consistent 
with nonergodicity and with the g-generalized Central Limit Theorem. It confirms the more-than- 
centennial prediction by J.W. Gibbs that standard statistical mechanics are not applicable for long- 
range interactions (i.e., for < a < 1) due to the divergence of the canonical partition function. 
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More than one century ago, in his historical book El- 
ementary Principles in Statistical Mechanics [1], J. W. 
Gibbs remarked that systems involving long-range inter- 
actions will be intractable within his and Boltzmann the- 
ory, due to the divergence of the partition function. This 
is of course the reason why no standard temperature- 
dependent thermostatistical quantities (e.g., a specific 
heat) can possibly be calculated for the free hydrogen 
atom, for instance. Indeed, unless a box surrounds the 
atom, an infinite number of excited energy levels accu- 
mulate at the ionization value, which yields a divergent 
canonical partition function at any finite temperature. 

To transparently extract the deep consequences of 
Gibbs' remark, in the present paper we focus on the in- 
fluence of the range of the interactions within an illustra- 
tive isolated classical system, namely the a-XY model [2] , 
whose Hamiltonian is given by 
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where the planar rotators are located at the sites of a 
d-dimcnsional hypercubic lattice with periodic boundary 
conditions. For d = 1, takes the values 1, 2, 3...; for 
d = 2, it takes the values 1, a/2, 2, . . . ; for d = 3, it takes 

the values 1, a/2, a/3, 2, The distance between any 

two rotators is taken to be the minimal one given the 
periodic boundary conditions. Without loss of general- 
ity we have considered unit moment of inertia, and unit 
first-neighbor coupling constant; pi and 9i are canonical 



conjugate pairs. At the fundamental state, all rotators 
are parallel, say &i = 0, Vz, which corresponds to the 
ferromagnetically fully ordered case. At high enough en- 
ergies, the values of {0i\ are randomly distributed, which 
corresponds to the paramagnetic phase. In between, 
a second order phase transition occurs. ^The potential 
energy per particle varies with N like N = ^^=2 

This quantity can be approximated, for a/d < 00, by 

N l-c/d_ 1 



d J ± drr r~ 
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which in turn behaves, 



when N -> 00, like N 1 -**^/^ - a/d) if < a/d < 1, 
like In A if a/d = 1, and like l/(a/d- 1) if a/d > 1. In 
other words, the total potential energy is extensive (in 
the thermodynamical sense) for a/d > 1, and nonexten- 
sive otherwise. In order to accommodate to a common 
practice, we might re- write the Hamiltonian % as follows: 
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which can now be considered as "extensive" for all values 
of a/d. This corresponds in fact to a rescaling of time 
(hence of pi), as shown in [2]. Also, this rewriting takes 
into account the fact that, for all values of a/d, the ther- 
modynamic energies (internal, Hclmholtz, Gibbs) grow 
like NN, the entropy, volume, magnetization, number of 
particles, etc, grow like N (i.e., remain extensive for both 
regions above and below a/d = 1), and the temperature, 
pressure, external magnetic field, chemical potential, etc, 
must be scaled with N in order to have finite equations of 
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states [3]. The particular case a = recovers the HMF 
model ([4] and references therein); the a —> oo model 
corresponds to first-neighbor interactions (whose d = 1 
case has been analytically studied [5]). 

In addition to the above, it has already been shown 
that the special value a/d = 1 also emerges dynami- 
cally. Indeed, for N — > oo and energies corresponding to 
the paramagnetic region, the largest Lyapunov exponent 
of the many-body system remains finite and positive for 
a/d > 1, whereas gradually vanishes for < a/d < 1. It 
vanishes like N~ K 7 where K,(a/d) decreases from a posi- 
tive value (close to 1/3) to zero when a/d increases from 
zero to 1, and remains zero for a/d > 1. It is interesting 
to emphasize that k does not independently depend on 
(a, d), but only on the ratio a/d [2, 6, 7]. Consistently 
with the fact that, for all values of the energy per particle 
u in the paramagnetic region, the Lyapunov exponents 
vanish in the limit N — > oo, k docs not depend on u. 

Let us briefly mention at this point that the breakdown 
of ergodicity which emerges for a/d < 1 [8] points to- 
wards the inadequacy of the Boltzmann-Gibbs (BG) the- 
ory. It is the aim of nonextensive statistical mechanics [3, 
9] to provide a way out of this kind of difficulty. Within 
this generalized theory, the stationary state is expected 
to yield a probability distribution p q = e q q /Z q ((3 q ) 
with Zqifiq) = J dpi...dpNd9i...d9Ne q ^ q ' H , where e q = 
[1 + (1 - q)x} 1 / { - 1 -^ (q G TZ; ef = e x ). The index q is 
expected to characterize universality classes, possibly a 
function q(a/d) to be different from 1 for < a/d < 1, 
and equal to 1 for a/d > 1. If this is so, an interesting 
quantity would of course be the one-momentum marginal 
probability P(pi) — J dpi-.-dpNdQ\...dQ^e q ^ q ' H ' /Z q . The 
functional form of P(pi) is unknown. A possibility could 
however be that, in the TV — > oo limit, we simply have 

P(pi) oc e ?r f ,mPl ^ 2 , i.e., a g m -Gaussian form, where m 
stands for momentum. Indeed, g-Gaussians emerge ex- 
tremely frequently in nonextensive- like systems (see, e.g., 
[10-15]; see also [16]). The index q m could depend not 
only on a/d, but also, in principle, on u (we remind 
that the d = 1 critical point for < a < 1 is known 
to be u c = 3/4). Now that we have outlined a possible 
thcrmostatistical scenario, let us present the molecular- 
dynamical results that we have obtained for the d = 1 
Hamiltonian (2) with fixed (N,u), the total energy be- 
ing Nu. We have used the Yoshida 4t/i-order symplcctic 
algorithm, and we have checked also through the stan- 
dard 4i/i-order Rungc-Kutta one. We present in Fig. 1 
the "temperature" T(t) = 2K(t)/N, where K{t) is the 
time-dependent total kinetic energy of Hamiltonian T-L. 
The class of initial conditions that we run are the so- 
called water-bag for the moments, with 0i = 0, Vi. As 
verified many times in the literature, a quasi-stationary 
state (QSS) exists for < a/d < 1 and u ~ 0.69, af- 
ter which a crossover is observed to a state whose tem- 
perature coincides with that analytically obtained within 
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FIG. 1: Time dependence of T(t) = 2K(t)/N for a water 
bag typical single initial condition for (it, iV) = (0.69, 200000) 
and typical values of a. The upper horizontal line, at Too = 
0.7114 . . ., corresponds to the (a, iV) —¥ (oo, oo) model at u = 
0.69 [5]. The middle (lower) horizontal line, at T ~ 0.475 
(T = 0.380), indicates the BG thermal equilibrium tempera- 
ture (the QSS base temperature, corresponding to zero mag- 
netization), at it = 0.69 and < a < 1. 



BG statistical mechanics [4]. The lifetime of this QSS 
appears to diverge with diverging N. It has been long 
thought that, after this crossover, the system consistently 
adopts a BG distribution in Gibbs V space, and there- 
fore a Maxwellian distribution for P{pi). The facts that 
we now present reveal a much more complex situation, 
where robust g„-Gaussians (or distributions numerically 
very close to them) emerge before the crossover (just be- 
fore for most realizations of the initial conditions, but 
also quite before for not few of them) and remain so for 
huge times (practically for ever); n stands for numerical. 
This unexpected phenomenon occurs for u both below 
and above u c = 3/4, and for a both below and above 
a = 1 (up to a ~ 2). Let us emphasize that these q n - 
Gaussians only develop their full shape if sufficient time 
has been run in order that the apparently stationary state 
has practically been attained. This time is extremely long 
for < u <C 3/4 because the system is then almost inte- 
grable (indeed, the Hamiltonian can be straightforwardly 
checked to become very close to that of N coupled har- 
monic oscillators, by using cos^ — 0j) ~ 1 — h (6i — 9j ) 2 ) , 
and is also extremely long for u ^ 3/4 because once again 
the system is almost integrable (indeed, the Hamiltonian 
can be straightforwardly checked to become now very 
close to N independent localized rotators). Let us de- 
tail now how the single-initial-condition one-momentum 
distributions P(p) are calculated within large time re- 
gions where T is sensibly constant: for each value of i, 
we register its pi at very many (noted n) successive times 
separated by an interval r, and then, following the recipe 
of the g-generalized Central Limit Theorem [17], calcu- 
late its arithmetic average pi (thus corresponding to the 
interval t G [t m i n ,t ma x] with t max - t min = nr). We 




FIG. 2: A typical single-initial-condition one-momentum 
distribution P(p) for N = 10 6 , u = 0.69, a = 0.9, r = 1 
(corresponding to 5 molecular-dynamical algorithmic steps), 
calculated in the region (tmin,t ma x) = (200000, 500000) (see 
Inset), where the temperature coincides with that analyt- 
ically calculated within BG statistical mechanics, namely 
T(oo) = 2K(oo)/N ~ 0.475. The total energy Nu is con- 
served within a relative precision of 10" 5 or better. The con- 
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tinuous curve corresponds to P(p)/Po = e q 

o(P ) 



with 

(?n, Pq„ a> ) = (1.58, 11.2). The value of q n corresponds to the 
red open circle in Fig. 3. Notice that here 1//3™ / T. Each 
distribution has been rescaled with its own Pq. 



then plot the histogram for the TV arithmetic averages, 
as illustrated in Fig. 2. All the histograms that we have 
obtained for sufficiently large times t are well fitted with 

e qf q " P ^ \ w hh (q n ,/3q n ) depending on (a, u, N, t) as well 
as on {t m i n , t max ). To check the quality of the fit we in- 
troduce (see Fig. 3) a conveniently q-generalized kurtosis 
(referred to as q-kurtosis), defined as follows: 



Ho dppHPip)} 2 "- 1 / JZv dp [P(.P)} 
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where we have used the escort distributions (see [18] and 
references therein). These distributions have the remark- 
able advantage of being finite up to q = 3, which is pre- 
cisely the value below which q-Gaussians are normaliz- 
able, i.e. dp P o eq 0qp2/2 = 1 (q < 3). The use of the 
standard kurtosis k\ = (p 4 )/3(p 2 ) 2 has the considerable 
disadvantage that (p 2 ) diverges for q > 5/3, and (p 4 ) di- 
verges for q > 3/2. Hence ki becomes useless for q > 3/2, 
and it happens that some of the distributions that we ob- 
serve do exhibit q n > 3/2. If we use a q- Gaussian P(p) 
within Eq. (3), we obtain, through a relatively easy cal- 
culation, K q (q) = fx 2 , as also shown in Fig. 3. 

In Figs. 4 and 5 we illustrate (q n ,l3q n ) as functions 
of (a,u,r) for large values of N. All the results for q n 
have been also reported in Fig. 3. One of the interesting 
features that we can observe is that in all cases q n ap- 
proaches the BG value q — 1 when r increases. However, 



FIG. 3: q n and g-kurtosis n q „ that have been obtained from 
the histograms corresponding to typical values of a (numbers 
indicated on top of the points). The red circle corresponds 
to Fig. 2. The continuous curve K q = (3 — g)/(l + q) is 
the analytical one obtained with g-Gaussians. Notice that 
k 9 is finite up to q = 3 (maximal admissible value for a q- 
Gaussian to be normalizable) , and that it does not depend 
on j3 qn . The visible departure from the dotted line at n q = 1 
corresponding to a Maxwellian distribution, neatly reflects 
the departure from BG thermostatistics. 



this approach is nearly exponential for (a < 1, u < 0.75), 
(a > 1, u > 0.75), and (a > 1, u < 0.75), whereas it 
is extremely slow for (a < 1, u > 0.75) (notice that, 
in the latter exhibits a zero slope with re- 

gard to t at t = 1), precisely the region where the 
largest Lyapunov exponent approaches zero with increas- 
ing N. This suggests the following nonuniform conver- 
gence: limjv-voo lim T _j. 00 q n (a, u, N, r) = 1 (Va), whereas 
linv-Kx, lini/v-Kx) qn (a, u, N, r) > 1 (for < a < 1). 
Lack of computational strength has not allowed us to di- 
rectly verify this conjecture. This leaves as an interesting 
open question whether linv^oo limAr^oo q n (a, u, N, r) re- 
covers limjv->oo qm{ot, u, N), where the latter would cor- 
respond to successive approximations for increasingly 
large TV. 

Summarizing, it has been observed for at least one 
decade that, for < a < 1, the longstanding QSSs of the 
present model exhibit anomalous distributions (Vlasov- 
like for some classes of initial conditions, and different, 
including g-Gaussian-shaped, ones for other classes) for 
the momenta of the rotators, whereas nothing particu- 
larly astonishing was expected to occur once the system 
had done the crossover to the (presumably stationary) 
state whose temperature coincides with that analytically 
obtained within the BG theory. The present results (ob- 
tained from first principles, i.e., using essentially noth- 
ing but Newton's law) neatly show that, if time is large 
enough so that the crossover has occurred (as illustrated 
in the Inset of Fig. 2), the situation is far more com- 
plex. Indeed, robust and longstanding g-Gaussian dis- 
tributions are numerically observed under a wide variety 
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FIG. 4: a-dependences of (qn,Pq n ) for (u, r, TV) — 
(0.69, 1.0, iV), where TV = 200000 (TV = 30000) for a > 0.6 
(a < 0.5), with n never smaller than 300000. We have ver- 
ified the existence of finite-size effects, in particular, for a 
above and close to unity, q n slowly decreases with increasing 
TV. Notice that T a (oo) = K(oo)/2N ~ 0.475 up to a ~ 1.35, 
where it starts increasing (red full circles), and, for a>l, ap- 
proaches the analytical value = 0.7114 . . .[5] (by using the 
values that we have obtained up to a = 40, we observe that 
approximatively — T a ~ 0.4/a 2 for a S> 1). The red open 
circles correspond to the example in Fig. 2 (also indicated in 
Fig. 3). The dependence of T on (a,t) is noted T a (t), hence 
confusion between T(oo) and Too must be avoided. The full 
(open) triangles have been obtained from rescaled histograms 
where the momenta have been divided by the standard devi- 
ation a (multiplied by Po, as illustrated in Fig. 2). The error 
bars corresponding to the triangles are of the same order; the 
error bars of T(oo) are of the order of the full circles (red). 
Naturally, Po x a is nearly constant; to take into account 
the numerical deviations (from a strict constant) due to pa- 
rameters such as (TV, n, r), we have normalized both fi q ^} and 
fiqn i n such a way that the analytical value T x — 0.7114 . . . 
is recovered. 
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FIG. 5: r-dependence of q n for TV = 200000, (train, t ma x) = 
(90000, 1090000) (hence n = 1000000 for r = 1), and typical 
values of u above and below the critical value u c = 0.75, and 
of a above and below the special value a = 1 (see [2]). All the 
errors bars are of the same order of those indicated on the red 
empty triangles. Inset: r-dependence of [g„ (r) — l]/[g n (l) — l]. 



of situations. The fact that the temperature be the one 
predicted within the BG theory appears to be necessary 
but not sufficient for standard statistical mechanics to be 
applicable [19]. Indeed, the shape of the momenta distri- 
butions can considerably differ from Gaussians, and it is 
only when the correlations become negligible (i.e., when 
t ^> 1 and/or a 1) that the classical Maxwellian distri- 
bution (with = T) is to be (numerically) recovered. 
This example shows the great thermostatistical richness 
that a breakdown of ergodicity can cause. It also serves as 
an invitation for deeper analysis of the thermal statistics 
of all those very many models in the literature that are 
definitively nonergodic (e.g., glasses, spin-glasses, among 
others), and for which, nevertheless, the BG theory is 
used without further justification. 
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